Qualitative Risk Analysis _ Inha
Kyungtak Kim
2020 3 26
# 패키지 설치
#install.packages("sf")
#install.packages("tmap")
#install.packages("dplyr")library(tidyverse)
library(sf)
library(tmap)
Sys.setenv(Language="En")
library(caret)
library(knitr)
library(leaflet)
library(rgdal)
library(htmltools)RISK = f(Hazard, Exposure, Vulnerability, Capacity)
# 원본 데이터 읽기
hazard <- read.csv('output/hazard_result.csv')
exposure <- read.csv('output/exposure_result.csv')
vulnerability <- read.csv('output/vulnerability_result.csv')
capacity <- read.csv('output/capacity_result.csv')
#데이터 결합
DB <- cbind(hazard, c(exposure[,4:5],vulnerability[,4:5],capacity[,4:6]))
head(DB, 3)## Name NameK SGG X16_hazard X17_hazard X18_hazard
## 1 Gangwon Gangneung-si 강원도 강릉시 42150 1.0000000 1.0000000 1.0000000
## 2 Gangwon Goseong-gun 강원도 고성군 42820 0.5264574 0.5463284 0.6153013
## 3 Gangwon Donghae-si 강원도 동해시 42170 0.6123409 0.6041215 0.5936511
## X16_exposure X17_exposure X16_vulnerability X17_vulnerability X16_capacity
## 1 0.3038772 0.3239064 0.5435461 0.4354101 0.7983780
## 2 0.1938716 0.1884191 0.6601949 0.6438191 0.7478007
## 3 0.3807366 0.3786649 0.3719887 0.3689875 0.7751020
## X17_capacity X18_capacity
## 1 0.8558765 0.8870795
## 2 0.8281501 0.9009704
## 3 0.8444816 0.8442552RISK 계산
# 16년~17년 홍수피해위험지수 산정
result_index_16 <- as.data.frame((rowSums(DB[,c("X16_hazard","X16_exposure","X16_vulnerability","X16_capacity")]))/4)
colnames(result_index_16) <- c("X16_result_index")
result_index_17 <- as.data.frame((rowSums(DB[,c("X17_hazard","X17_exposure","X17_vulnerability","X17_capacity")]))/4)
colnames(result_index_17) <- c("X17_result_index")
result_index <- cbind(DB[,1:3], c(result_index_16,result_index_17))
# 홍수피해위험지수 표준화 함수 설정
standard <- function(x){
return((x-min(x))/(max(x)-min(x)))
}
# 연도별 데이터 프레임에 표준화 적용
result <- as.data.frame(lapply(result_index[,4:5],standard))
colnames(result) <- c("X16_result", "X17_result")
result <- cbind(DB[,1:3], result)
head(result, 3)## Name NameK SGG X16_result X17_result
## 1 Gangwon Gangneung-si 강원도 강릉시 42150 0.8100104 0.9003977
## 2 Gangwon Goseong-gun 강원도 고성군 42820 0.4881184 0.6160862
## 3 Gangwon Donghae-si 강원도 동해시 42170 0.4954857 0.6088050연도별 확률밀도함수
result_p <- result %>%
select(-Name)
head(result_p, 3)## NameK SGG X16_result X17_result
## 1 강원도 강릉시 42150 0.8100104 0.9003977
## 2 강원도 고성군 42820 0.4881184 0.6160862
## 3 강원도 동해시 42170 0.4954857 0.6088050result_p_p <- result_p %>% # pivoting
pivot_longer(c("X16_result", "X17_result"),
names_to = "year",
values_to = "result")
result_p_p %>%
ggplot()+
geom_density(aes(x=result, y=..density.., color=year))
result_p %>%
ggplot(aes(X17_result))+
geom_histogram(bins=100)
result_p_p %>%
group_by(NameK) %>%
mutate(mean=mean(result))%>%
ggplot(aes(x=fct_reorder(NameK, mean),
y=result))+
geom_boxplot()+
coord_flip()
result_p_p %>%
group_by(NameK) %>%
mutate(mean=mean(result))%>%
filter(mean < 0.35) %>% #35% 이하
ggplot(aes(x=fct_reorder(NameK, mean),
y=result))+
geom_boxplot()+
coord_flip()
result_p_p %>%
group_by(NameK) %>%
mutate(mean=mean(result))%>%
filter(mean > 0.75) %>% #75% 이상
ggplot(aes(x=fct_reorder(NameK, mean),
y=result))+
geom_boxplot()+
coord_flip()
result_p %>%
mutate(dif=(X17_result - X16_result)) %>%
filter(NameK == "서울특별시")## NameK SGG X16_result X17_result dif
## 1 서울특별시 11000 0.7166389 0.7321177 0.01547886result_p_dif <- result_p%>%
mutate(dif=(X17_result - X16_result)) %>%
arrange(-dif)
knitr::kable(result_p_dif[1:10, ]) # 침수구역내 총인구가 늘어난 시군| NameK | SGG | X16_result | X17_result | dif |
|---|---|---|---|---|
| 강원도 양양군 | 42830 | 0.6178257 | 0.7966093 | 0.1787836 |
| 경상남도 함안군 | 48730 | 0.4164833 | 0.5869667 | 0.1704833 |
| 강원도 평창군 | 42760 | 0.7611326 | 0.9193971 | 0.1582644 |
| 강원도 춘천시 | 42110 | 0.5502588 | 0.7025579 | 0.1522991 |
| 경상북도 영양군 | 47760 | 0.5093570 | 0.6605020 | 0.1511450 |
| 경상북도 군위군 | 47720 | 0.4829632 | 0.6289233 | 0.1459600 |
| 강원도 영월군 | 42750 | 0.4185902 | 0.5598876 | 0.1412974 |
| 경상북도 의성군 | 47730 | 0.4869166 | 0.6280056 | 0.1410890 |
| 강원도 고성군 | 42820 | 0.4881184 | 0.6160862 | 0.1279678 |
| 강원도 횡성군 | 42730 | 0.3402637 | 0.4658988 | 0.1256351 |
knitr::kable(result_p_dif[152:161, ]) # 침수구역내 총인구가 줄어든 시군| NameK | SGG | X16_result | X17_result | dif | |
|---|---|---|---|---|---|
| 152 | 경기도 안산시 | 41270 | 0.4316604 | 0.4016536 | -0.0300069 |
| 153 | 충청남도 청양군 | 44790 | 0.7362989 | 0.7017132 | -0.0345857 |
| 154 | 경기도 연천군 | 41800 | 0.6509595 | 0.6125975 | -0.0383620 |
| 155 | 경기도 부천시 | 41190 | 0.4278185 | 0.3781803 | -0.0496382 |
| 156 | 충청남도 부여군 | 44760 | 0.7773486 | 0.7275116 | -0.0498370 |
| 157 | 인천광역시 | 28000 | 0.5263018 | 0.4756885 | -0.0506133 |
| 158 | 경기도 시흥시 | 41390 | 0.2061888 | 0.1555098 | -0.0506790 |
| 159 | 충청북도 청주시 | 43110 | 0.3649189 | 0.3043571 | -0.0605618 |
| 160 | 경기도 고양시 | 41280 | 0.7098045 | 0.6421808 | -0.0676237 |
| 161 | 충청북도 증평군 | 43745 | 0.3408069 | 0.2562014 | -0.0846056 |
result_p_p %>%
group_by(year) %>%
ggplot(aes(result, SGG))+
geom_point(aes(color=factor(SGG)))+
facet_grid(. ~year)+
theme(legend.position = "none")
Mapping
# 시군 shp 파일 불러오기
analysis <- st_read("input/analysis.shp")## Reading layer `analysis' from data source `C:\00_R\0_Git\KRM_inha\input\analysis.shp' using driver `ESRI Shapefile'
## Simple feature collection with 161 features and 3 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 746109.3 ymin: 1458771 xmax: 1387956 ymax: 2068444
## proj4string: +proj=tmerc +lat_0=38 +lon_0=127.5 +k=0.9996 +x_0=1000000 +y_0=2000000 +ellps=GRS80 +units=m +no_defs# 폴리곤 에러 체크(기본 파일을 에러 수정한 파일로 변경하였음)
#st_is_valid(analysis)
#library(lwgeom)
#analysis <- st_make_valid(analysis)
st_is_valid(analysis)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [136] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [151] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE# shp파일에 연도별 홍수피해위험지수(표준화 적용) 추가
analysis <- right_join(analysis, result[,3:5])## Joining, by = "SGG"# 폴리곤 단순화
analysis_simp <- st_simplify(analysis, dTolerance = 50)
# 결과 확인
tmap_mode("plot")## tmap mode set to plottingbreaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)
facets=c("X16_result", "X17_result")
tm_shape(analysis_simp)+
tm_polygons(facets,
breaks=breaks,
palette = c("green", "greenyellow", "yellow", "orange", "red"),
legend.reverse = TRUE)+
tm_facets(ncol = 2)+
tm_layout(legend.position = c("right", "bottom"))+
tm_compass(type = "rose",
position = c("right", "top"),
size = 2.0)+
tm_scale_bar(breaks = c(0, 25, 50, 100, 150, 200),
position = c("left", "bottom"))
######################
#library(leaflet)
#library(rgdal)
#library(htmltools)a <- st_transform(analysis_simp, 4326)
pal <- colorBin(
palette=c("green", "greenyellow", "yellow", "orange", "red"),
domain=NULL,
bins = c(0, .2, .4, .6, 0.8, 1),
pretty = FALSE)
leaflet(a) %>%
setView(lng = 128, lat = 35.9, zoom = 7) %>%
# base groups
addPolygons(color = ~pal(X16_result),
weight = 1,
smoothFactor = 0.5,
opacity = 1.0,
fillOpacity = 0.5,
label = ~htmlEscape(NameK),
popup = ~htmlEscape(X16_result),
highlightOptions = highlightOptions(color = "white",
weight = 2,
bringToFront = TRUE),
group="result 2016") %>%
addPolygons(color = ~pal(X17_result),
weight = 1,
smoothFactor = 0.5,
opacity = 1.0,
fillOpacity = 0.5,
label = ~htmlEscape(NameK),
popup = ~htmlEscape(X17_result),
highlightOptions = highlightOptions(color = "white",
weight = 2,
bringToFront = TRUE),
group="result 2017") %>%
#overlay groups
addProviderTiles(providers$Esri.WorldStreetMap,
group="Esri") %>%
addProviderTiles(providers$CartoDB.Positron,
group="CartoDB") %>%
addLegend("bottomright",
pal = pal,
values = ~X17_result,
title = "RISK Index",
labFormat = labelFormat(digits=10),
opacity = 1) %>%
hideGroup("CartoDB") %>%
#Layer controls
addLayersControl(baseGroups = c("result 2016", "result 2017"),
overlayGroups = c("Esri", "CartoDB"),
options=layersControlOptions(collapsed=FALSE))#############################결과값 저장
write.csv(result, 'output/final_result.csv')
# 열 명칭별 의미
# Name : 161개 시군별 영문명
# NameK : 161개 시군별 한글명
# SGG : 시군구 코드
# X16_hazard : 16년도 hazard 지수(표준화 적용)
# X17_hazard : 17년도 hazard 지수(표준화 적용)
# X18_hazard : 18년도 hazard 지수(표준화 적용)
# X16_exposure : 16년도 Exposure 지수(표준화 적용)
# X17_exposure : 17년도 Exposure 지수(표준화 적용)
# X16_vulnerability : 16년도 Vulnerability 지수(표준화 적용)
# X17_vulnerability : 17년도 Vulnerability 지수(표준화 적용)
# X16_capacity : 16년도 Capacity 지수(표준화 적용)
# X17_capacity : 17년도 Capacity 지수(표준화 적용)
# X18_capacity : 18년도 Capacity 지수(표준화 적용)
# X16_result_index : 16년도 홍수피해위험지수
# X17_result_index : 17년도 홍수피해위험지수
# X16_result : 16년도 홍수피해위험지수(표준화 적용)
# X17_result : 17년도 홍수피해위험지수(표준화 적용)